##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
library(RColorBrewer)
library(GSVA)

##########################################################################################

option_list <- list(
    make_option(c("--ddr_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/RSEM/CombineTPM.tsv"
    ddr_file <- "~/WGS_Scripts/v2/PathWay/DDR.tsv"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/ddr_gsva"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

ddr_file <- opt$ddr_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file

dir.create(out_path , recursive = T)

##########################################################################################

dat_tpm <- data.frame(fread(rsem_file))
colnames(dat_tpm) <- gsub( "X" , "" , colnames(dat_tpm))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

dat_tpm$gene_id <- sapply( strsplit(dat_tpm$gene_id , "[.]") , "[" , 1 )
dat_tpm_all_out <- merge( dat_gtf , dat_tpm)

#去除重名基因
dat_tpm_all_out <- dat_tpm_all_out[!duplicated(dat_tpm_all_out$Hugo_Symbol),]
dat_tpm_all_out <- dat_tpm_all_out[,-1]
colnames(dat_tpm_all_out)[1] <- "gene_id"
rownames(dat_tpm_all_out) <- dat_tpm_all_out$gene_id
dat_tpm_all_out <- dat_tpm_all_out[,-1]
dat_tpm_all_out <- as.matrix(dat_tpm_all_out)

##########################################################################################
## log转化
dat <- as.matrix(log2(dat_tpm_all_out+1))

##########################################################################################

dat_pathway <- fread(ddr_file)

##########################################################################################

pathlist <- unique(sapply(strsplit(dat_pathway$Pathway , ";"), "[" , 1))
all_path <- list()
for( path in pathlist ){
    geneName <- dat_pathway$GeneName[grep(path,dat_pathway$Pathway)]
    all_path[[path]] <- geneName
}

all_path[["All_DDR"]] <- unique(dat_pathway$GeneName)

##########################################################################################
## "Gaussian" for logCPM,logRPKM,logTPM, "Poisson" for counts
gsva.es <- gsva(dat, all_path , kcdf = "Gaussian" , verbose = FALSE)

## 标化到0-1
#gsva_out <- data.frame(apply( gsva.es , 1 , function(x) {x=(x-min(x))/(max(x)-min(x))} ))
gsva_out <- data.frame(apply( gsva.es , 1 , function(x) {x} ))
gsva_out$Sample <- rownames( gsva_out )
gsva_out <- gsva_out[,c(ncol(gsva_out) , 1:(ncol(gsva_out)-1))]

out_name <- paste0(out_path , "/gsva_ddr.tsv")
write.table( gsva_out , out_name , row.names = F , sep = "\t" , quote = F )